/*
  several function used for PbPb combined spectra
  Blast Wave is also implemented here
  further documentation will come
  
  author: Roberto Preghenella
  email : preghenella@bo.infn.it
*/

//http://alisoft.cern.ch/viewvc/trunk/PWGLF/SPECTRA/PiKaPr/COMBINED/SpectraUtils.C?view=markup&root=AliRoot                                                                                                                                                  

const Double_t fxmax = 25;
Double_t fPtMax[]={100,100,100};
Double_t fPtMin[]={0,0,0};
Double_t fChi2 = -999;
Double_t fNdata = 0;

void
BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);

/*****************************************************************/
/* BOLTZMANN-GIBBS BLAST-WAVE */
/*****************************************************************/

static TF1 *fBGBlastWave_Integrand = NULL;
Double_t
BGBlastWave_Integrand(const Double_t *x, const Double_t *p)
{
  
  /* 
     x[0] -> r (radius)
     p[0] -> mT (transverse mass)
     p[1] -> pT (transverse momentum)
     p[2] -> beta_max (surface velocity)
     p[3] -> T (freezout temperature)
     p[4] -> n (velocity profile)
  */
  
  Double_t r = x[0];
  Double_t mt = p[0];
  Double_t pt = p[1];
  Double_t beta_max = p[2];
  Double_t temp_1 = 1. / p[3];
  Double_t n = p[4];

  Double_t beta = beta_max * TMath::Power(r, n);
  if (beta > 0.9999999999999999) beta = 0.9999999999999999;
  Double_t rho = TMath::ATanH(beta);
  Double_t argI0 = pt * TMath::SinH(rho) * temp_1;
  if (argI0 > 700.) argI0 = 700.;
  Double_t argK1 = mt * TMath::CosH(rho) * temp_1;
  //  if (argI0 > 100 || argI0 < -100)
  //    printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", r, pt, beta_max, 1. / temp_1, n, mt, beta, rho, argI0, argK1);
  return r * mt * TMath::BesselI0(argI0) * TMath::BesselK1(argK1);
  
}

Double_t
BGBlastWave_Func_OneOverPt(const Double_t *x, const Double_t *p)
{
  /* 1/pt dN/dpt */
  
  Double_t pt = x[0];
  Double_t mass = p[0];
  Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
  Double_t beta_max = p[1];
  Double_t temp = p[2];
  Double_t n = p[3];
  Double_t norm = p[4];
  
  if (!fBGBlastWave_Integrand)
    fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
  fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
  Double_t integral = fBGBlastWave_Integrand->Integral(0., 1.);

  //return norm * integral;
  //now it is not OneOverPt but just dN/dpt
  return norm * integral *pt;
}


TF1 * BGBlastWave_OneOverPT(const Char_t *name, Double_t mass, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm = 1.e6)
{
  
  TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func_OneOverPt, 0., fxmax, 5);
  fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
  fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
  fBGBlastWave->FixParameter(0, mass);
  //fBGBlastWave->SetParLimits(1, 0.01, 0.99);
  //fBGBlastWave->SetParLimits(2, 0.01, 1.);
  //fBGBlastWave->SetParLimits(3, 0.01, 10.);
  return fBGBlastWave;
}

/*****************************************************************/
/*****************************************************************/
/*****************************************************************/

Int_t nBW;
TF1 *fBGBW[1000];
TGraphAsymmErrors *gBW[1000];

TList *
BGBlastWave_GlobalFit(TString tag, TObjArray *data, Double_t *mass, const Double_t *tmpPtMin, const Double_t *tmpPtMax, TList *lout, Double_t profile = 1, Bool_t fixProfile = kFALSE)
{

  printf("Please check fixProfile %d profile %e\n", fixProfile, profile);

  /* get data */
  nBW = data->GetEntries();
  for (Int_t idata = 0; idata < nBW; idata++) {
    gBW[idata] = (TGraphAsymmErrors *)data->At(idata);
    gBW[idata]->SetName(Form("gBW%d", idata));
  }

  /* init BG blast-wave functions */
  for (Int_t idata = 0; idata < nBW; idata++) {
    fPtMin[idata]=tmpPtMin[idata];
    fPtMax[idata]=tmpPtMax[idata];
    printf("init BG-BlastWave function #%d: mass = %f ptmin = %f ptmax = %f\n", idata, mass[idata], fPtMin[idata],  fPtMax[idata]);
    printf("init BG-BlastWave function #%d: mass = %f\n", idata, mass[idata]);
    fBGBW[idata] = BGBlastWave_OneOverPT(Form("fBGBW%d", idata), mass[idata]);
  }

  /* display data */
  TCanvas *cBW = new TCanvas("cBW");
  cBW->Divide(nBW, 1);
  for (Int_t idata = 0; idata < nBW; idata++) {
    cBW->cd(idata + 1);
    gBW[idata]->Draw("ap*");
  }
  cBW->Update();

  /* init minuit: nBW normalizations + 3 (beta, T, n) BG-BlastWave params */
  const Int_t nbwpars = 3;
  const Int_t nfitpars = nBW + nbwpars;
  TMinuit *minuit = new TMinuit(nfitpars);
  minuit->SetFCN(BGBlastWave_FCN);
  Double_t arglist[10];
  Int_t ierflg = 0;
  arglist[0] = 1;
  minuit->mnexcm("SET ERR", arglist, 1, ierflg);
  for (Int_t idata = 0; idata < nBW; idata++)
    minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);

  //with limits
  minuit->mnparm(nBW + 0, "<beta>", 0.65, 0.01, 0., 1., ierflg);
  minuit->mnparm(nBW + 1, "T", 0.1, 0.01, 0., 1., ierflg);
  minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
 
  if (fixProfile) minuit->FixParameter(nBW + 2);

  const Double_t tol = 1;
  /* set strategy */
  arglist[0] = 1;
  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);

  /* start MIGRAD minimization */
  arglist[0] = 500000;
  arglist[1] = tol;
  //minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
  minuit->mnexcm("MINIMIZE", arglist, 2, ierflg);
  printf("STRA 1 ierr %d\n", ierflg);

  /* set strategy */
  arglist[0] = 2;
  minuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);

  /* start MIGRAD minimization */
  arglist[0] = 500000;
  arglist[1] = tol;
  //minuit->mnexcm("MIGRAD", arglist, 2, ierflg);
  minuit->mnexcm("MINIMIZE", arglist, 2, ierflg);
  printf("STRA 2 ierr %d fChi2 %e\n", ierflg, fChi2);

  const Double_t finalChi2 = fChi2;
  printf("%s result MINIMIZE fChi2 %e\n", tag.Data(), fChi2);

  /* start IMPROVE minimization */
  arglist[0] = 500000;
  minuit->mnexcm("IMPROVE", arglist, 1, ierflg);
  //ierr 4:  IMPROVE HAS RETURNED TO REGION OF ORIGINAL MINIMUM
  //also IMPROVE seems to make failing MIGRAD convergent again for STAR 2 ierr 4.
  printf("IMPROVE ierr %d fChi2 %e\n", ierflg, fChi2);
  printf("%s result IMPROVE fChi2 %e\n", tag.Data(),  fChi2);

  /* start MINOS */
  arglist[0] = 500000;
  arglist[1] = nBW + 1;
  arglist[2] = nBW + 2;
  arglist[3] = nBW + 3;
  minuit->mnexcm("MINOS", arglist, 4, ierflg);
  printf("MINOS ierr %d\n", ierflg);

  /* print results */
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
  minuit->mnprin(4, amin);

  /* get parameters */
  Double_t beta, betae, betaeplus, betaeminus, betagcc, temp, tempe, tempeplus, tempeminus, tempgcc, prof, profe, profeplus, profeminus, profgcc;
  minuit->GetParameter(nBW + 0, beta, betae);
  minuit->mnerrs(nBW + 0, betaeplus, betaeminus, betae, betagcc);
  minuit->GetParameter(nBW + 1, temp, tempe);
  minuit->mnerrs(nBW + 1, tempeplus, tempeminus, tempe, tempgcc);
  minuit->GetParameter(nBW + 2, prof, profe);
  minuit->mnerrs(nBW + 2, profeplus, profeminus, profe, profgcc);
  Double_t beta_max = 0.5 * (2. + prof) * beta;
  Double_t norm[1000], norme[1000];
  for (Int_t idata = 0; idata < nBW; idata++)
    minuit->GetParameter(idata, norm[idata], norme[idata]);

  const Double_t nfree = minuit->GetNumFreePars();
  const Double_t ndof = fNdata - nfree;
  const Double_t redChi2 = finalChi2/ndof;
  printf("%s result nfree %.0f fNdata %.0f ndof %.0f MINIMIZE finalChi2 %e redChi2 %e\n", tag.Data(), nfree, fNdata, ndof, finalChi2, redChi2);

  /* printout */
  printf("*********************************\n");
  printf("%s result beta_max = %f\n", tag.Data(), beta_max);
  printf("%s result <beta>   = %f +- %f (e+ = %f, e- = %f)\n", tag.Data(), beta, betae, betaeplus, betaeminus);
  printf("%s result T        = %f +- %f (e+ = %f, e- = %f)\n", tag.Data(), temp, tempe, tempeplus, tempeminus);
  printf("%s result n        = %f +- %f (e+ = %f, e- = %f)\n", tag.Data(), prof, profe, profeplus, profeminus);

  /* 1-sigma contour */
  minuit->SetErrorDef(1);
  TGraph *gCont1 = NULL;
  gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
  if (gCont1) gCont1->SetName("gCont1");
  printf("cont1 ierr %d\n", gCont1?0:4);

  /* 2-sigma contour */
  minuit->SetErrorDef(4);
  TGraph *gCont2 = NULL;
  //  gCont2 = (TGraphAsymm *) minuit->Contour(50, nBW + 0, nBW + 1);
  if (gCont2) gCont2->SetName("gCont2");

  /* display fits */
  printf("calculating gfit and gratio...\n");
  TGraphAsymmErrors * gfit[10]={0x0,0x0,0x0};
  TGraphAsymmErrors * gratio[10]={0x0,0x0,0x0};

  for (Int_t idata = 0; idata < nBW; idata++) {
    cBW->cd(idata + 1);

    fBGBW[idata]->SetParameter(4, norm[idata]);
    fBGBW[idata]->SetParameter(1, beta_max);
    fBGBW[idata]->SetParameter(2, temp);
    fBGBW[idata]->SetParameter(3, prof);
    //fBGBW[idata]->Draw("same");

    gfit[idata]=new TGraphAsymmErrors;
    gfit[idata]->SetName(Form("gfit%d", idata));

    gratio[idata]=new TGraphAsymmErrors;
    gratio[idata]->SetName(Form("gratio%d", idata));

    for(Int_t ip = 0; ip< gBW[idata]->GetN(); ip++){
      const Double_t xx = gBW[idata]->GetX()[ip];
      const Double_t fitY = fBGBW[idata]->Eval(xx);

      gfit[idata]->SetPoint(ip, xx, fitY);
      printf("setting point for gfit %d %d: %e %e\n", idata,  ip, xx, fitY);

      const Double_t EXlow = gBW[idata]->GetEXlow()[ip];
      const Double_t EXhigh = gBW[idata]->GetEXhigh()[ip];

      const Double_t dataY = gBW[idata]->GetY()[ip];
      const Double_t dataEYlow = gBW[idata]->GetEYlow()[ip];
      const Double_t dataEYhigh = gBW[idata]->GetEYhigh()[ip];

      if(fitY/dataY > 0.01){
        const Int_t nratio = gratio[idata]->GetN();
        gratio[idata]->SetPoint(nratio, xx, dataY/fitY);
        gratio[idata]->SetPointError(nratio, EXlow, EXhigh, dataEYlow/fitY, dataEYhigh/fitY);
        printf("setting point for gratio %d %d: %e %e %e %e\n", idata,  nratio, xx, dataY/fitY, dataEYlow/fitY, dataEYhigh/fitY);
      }
    }
    gfit[idata]->SetLineColor(kRed);
    gfit[idata]->Draw("l");
  }
  cBW->Update();

  /* histo params */
  printf("setting hBW...\n");
  TH1D *hBW = new TH1D("hBW", "", 3, 0., 3.);
  hBW->SetBinContent(1, beta);
  hBW->SetBinError(1, betae);
  hBW->SetBinContent(2, temp);
  hBW->SetBinError(2, tempe);
  hBW->SetBinContent(3, prof);
  hBW->SetBinError(3, profe);

  /* BW graph */
  printf("setting gBetaT...\n");
  TGraphAsymmErrors *gBetaT = new TGraphAsymmErrors();
  gBetaT->SetName("gBetaT");
  gBetaT->SetPoint(0, beta, temp);
  gBetaT->SetPointEXlow(0, TMath::Abs(betaeminus));
  gBetaT->SetPointEXhigh(0, TMath::Abs(betaeplus));
  gBetaT->SetPointEYlow(0, TMath::Abs(tempeminus));
  gBetaT->SetPointEYhigh(0, TMath::Abs(tempeplus));

  printf("\nputting to TList...\n\n");

  /* prepare output array */
  TList *outoa = new TList;
  outoa->SetName("outoa");
  for (Int_t idata = 0; idata < nBW; idata++) {
    outoa->Add(gBW[idata]);
    outoa->Add(gfit[idata]);
    outoa->Add(gratio[idata]);
    //outoa->Add(fBGBW[idata]);
  }
  outoa->Add(hBW);
  outoa->Add(gBetaT);
  if (gCont1) outoa->Add(gCont1);
  if (gCont2) outoa->Add(gCont2);
  outoa->Add(cBW);

  //for lout
  lout->Add(gBetaT);
  if (gCont1) lout->Add(gCont1);

  return outoa;
}

void 
BGBlastWave_FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{

  /* beta -> beta_max */
  Double_t beta = par[nBW+0];
  Double_t T = par[nBW+1];
  Double_t n = par[nBW+2];
  Double_t beta_max = 0.5 * (2. + n) * beta;
#if 0
  /* check beta_max */
  if (beta_max >= 1. || beta_max <= 0.) {
    f = kMaxInt;
    return;
  }
  /* check T */
  if (T <= 0.) {
    f = kMaxInt;
    return;
  }
#endif

  Double_t pt, /*pte,*/ val, vale, func, pull, chi = 0;
  /* loop over all the data */
  fNdata = 0;
  for (Int_t iBW = 0; iBW < nBW; iBW++) {
    /* set BGBW parameters */
    fBGBW[iBW]->SetParameter(4, par[iBW]);
    fBGBW[iBW]->SetParameter(1, beta_max);
    fBGBW[iBW]->SetParameter(2, T);
    fBGBW[iBW]->SetParameter(3, n);
    /* loop over all the points */
    for (Int_t ipt = 0; ipt < gBW[iBW]->GetN(); ipt++) {
      pt = gBW[iBW]->GetX()[ipt];
      //introduce a pt cut off
      if(pt>fPtMax[iBW] || pt<fPtMin[iBW]){
        continue;
      }

      fNdata++;

      //pte = gBW[iBW]->GetEX()[ipt];
      val = gBW[iBW]->GetY()[ipt];
      vale = gBW[iBW]->GetEYlow()[ipt];
      func = fBGBW[iBW]->Eval(pt);
      //      func = fBGBW[iBW]->Integral(pt - pte, pt + pte);
      pull = (val - func) / vale;
      chi += pull * pull;
    }
  }

  f = chi;

  fChi2 = f;
}

